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Abstract 

Single nucleotide polymorphisms (SNPs) have become the marker of choice for genome-wide association studies in many 
species. High-throughput sequencing of RNA was developed primarily to analyze global gene expression, while it is an 
efficient way to discover SNPs from the expressed genes. In this study, we conducted transcriptome sequencing of the 
swimbladder of Takifugu rubripes using lllumina HiSeq2000 platform to identify gene-associated SNPs in the swimbladder. A 
total of 30,312,181 unique-mapped-reads were obtained from 44,736,850 raw reads. A total of 62,270 putative SNPs were 
discovered, which were located in 1 1,306 expressed genes and 2,246 scaffolds. The average minor allele frequency (MAP) of 
the SNPs was 0.26. GO and KEGG pathway analysis were conducted to analyze the genes containing SNPs. Validation of 
selected SNPs revealed that 54% of SNPs (26/48) were true SNPs. The results suggest that RNA-Seq is an efficient and cost- 
effective approach to discover gene-associated SNPs. In this study, a large number of SNPs were identified and these data 
will be useful resources for population genetic study, evolution analysis, resource assessment, genetic linkage analysis and 
genome-wide association studies. 
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Introduction 

Next-generation sequencing-based RNA-Seq analyses have 
dramatically changed the way to investigate the functional 
complexity of transcriptome in many organisms [1,2]. RNA-Seq 
approach is powerful for unraveling transcriptome complexity, 
identification of genes, gene-associated markers, regulatory non- 
coding RNAs and for alternative splicing analysis and expression 
profiling [3-5]. Transcriptome analysis using the next generation 
sequencing technologies have been widely reported in many 
species, including several aquaculture species such as catfish [6-8] , 
Atlantic cod [9], silver carp [10], pearl oyster [11], carp [12], and 
Amur ide [13]. 

Recently, RNA-Seq has also been used as an efficient and cost- 
effective method to comprehensively identify SNPs from tran- 
scribed regions in the genomes of several fish species. By 
sequencing of the pooled RNA samples from multiple individuals 
of channel catfish and blue catfish, a set of quality SNPs were 
identified including 342,104 intra-specific SNPs for channel 
catfish, 366,269 intra-specific SNPs for blue catfish, and 420,727 
inter-specific SNPs between channel catfish and blue catfish [6]. 
Similarly in carp, a total of 712,042 intra-stain SNPs were 
discovered in four strains, including mirror carp (483,276 SNPs), 
purse red carp (486, 629SNPs), Xingguo red carp (478,028 SNPs), 
and Yellow River carp (488,281 SNPs) [14]. Large sets of SNPs 
have also been reported in some other aquaculture species, such as 



the Eastern oyster [15], Atiantic salmon [16], Adantic cod [9] and 
rainbow trout [1 7] . 

Takifugu rubripes, widely distributed in the Asia, is one of the most 
important aquaculture species in China. In our laboratory, some 
SNPs makers associated with growth traits have been identified 
from the growth-related genes including Leptin, Mekmocortin 4 
Receptor {MC4R), Insulin-like growth factor {IGF), Myogenic factor 5 
{MjfSj, Growth hormone releasing hormone [GHRH], Myogenic factor 6 
{Myf6) [18]. Other genetic and genomic studies were also 
conducted with the focus on identification and characterization 
of microsatellite markers [19,20], construction of bacterial artificial 
chromosome (BAG) and expressed sequence tag (EST) library 
[21]. In addition to its importance in aquaculture, T. rubripes h also 
widely used as a model system in many scientific fields, especially 
in the evolutionary studies. The fugu genome has been completed, 
which is among the smallest vertebrate genomes. It has proven to 
be a useful Reference' genome for identifying genes and other 
functional elements in human and other vertebrate genomes, and 
for understanding the structure and evolution of vertebrate 
genomes [22-24]. 

The swimbladder in teleost fish is a specialized organ that 
regulates buoyancy. The homology of the fish swimbladder and 
mammalian lung has been well recognized based on morpholog- 
ical and embryological evidence. However, the molecular evidence 
of homology of swimbladder and the mammalian lung was not 
sufficient [25-27]. A large set of SNPs from the swimbladder 
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transcriptome of T. mhripe.s should provide valuable resources for 
swimbladder research, lung research and evolution studies of fish 
swimbladder and mammalian lung. 

In this study, we sequenced the transcriptome of the 
swimbladder of T. rubripes using lUumina HisSeq2000 platform 
to identify gene-associated SNPs. A total of 62,270 putative SNPs 
were discovered, which were located in 11,430 genes and 1,612 
scaffolds, and the average minor allele frequency (MAP) was 0.26. 
These SNPs should provide useful resources for evolution, 
population genetic study, resource assessment, genetic linkage 
analysis and genome-wide association studies. 

Results and Discussion 

Transcriptome sequencing 

lUumina sequencing was conducted to generate short sequence 
reads from the swimbladder of T. rubripes. A total of 30,312,181 
unique-mapped-reads were obtained from 44,736,850 raw reads 
after being mapped to the fugu T. rubripes fifth genome assembly 
from Ensembl database. The genome distribution of the uniquely 
mapped reads was assessed based on the RefSeq-defined gene 
models. As expected, the majority of reads (60%) were mapped 
onto exonic regions, while a large propotion of reads were mapped 
onto intergenic regions (Table 1). Similar observations have been 
reported in the studies of mouse and Caenorhabditis elegans [28,29] . 
The RNA-Seq data in this study has been deposited in the NCBI 
SRA database with the accession number of SRR 102267 7. 

SNP identification 

Compared with the fugu genome, a total of 62,270 putative 
SNPs were identified. The detailed SNP information was provided 
in Table SI. Of which, the number of homozygotes was 9,518 and 
the number of heterozygotes was 52,752. In these heterozygotes, 
the C/T and A/G were the most common types. In contrast, G/ 
T, A/C, G/C and A/T were the lease common types (Table 2). 

The SNPs were classified into several categories based on their 
locations in the genome, including inter-genic, down_stream (+lk), 
exon, intron, and up_stream (—1 k). As shown in Table 3, of the 
62,270 putative SNPs, 24,525 SNPs (39.38%) were identified in 
exons, which were highly represented, while 4,210 SNPs (6.76%) 
were identified in the introns, which were lowly represented. 

Minor allele frequency distribution 

Minor allele frequency (MAF) is an important factor for SNP 
locus evaluation. MAFs of SNPs were calculated from the 
sequence data. As shown in Figure 1 , the majority of SNPs have 
sequence derived minor allele frequencies ranging from 21% to 
25%, and the average MAF was 26% in putative SNPs identified 
from the swimbladder of T. rubripes. 



Table 1. The genome distribution of the mapped reads. 







Read distribution 


Number of reads 


Percentage 


Exonic region 


18,120,867 


59.78% 


Intergenic region 


9,776,865 


32.25% 


Intronic region 


1,435,258 


4.73% 


Exon-intron junction 


176,167 


0.58% 



doi:l 0.1 371 /journal.pone.0092502.t001 



Table 2. Summary of SNP types identified from the J. rubirpes 
swimbladder. 




SNP type 


Number 


Homozygote 


A 


1,887 


c 


2,914 


G 


2,887 


T 


1,830 


Heterozygote 


G/T 


4,730 


A/C 


4,725 


A/G 


16,972 


G/C 


4,665 


A/T 


4,586 


C/T 


17,074 


Total 


62,270 


doi:l 0.1 371 /journal.pone.0092502.t002 


Table 3. Classification of putative SNPs. 




SNP classification 


Number of putative SNPs 


Inter-genic 


12,903 


Down_stream{+l k) 


12,303 


Exon 


24,525 


Intron 


4,210 


Up_stream( — 1 k) 


8,329 


Total 


62,270 


Inter-genic SNPs were identified from regions between genes, while 
Down_stream{-t-l k) and Up_stream( — 1 k} represents SNPs identified from 
regions of 1 kb downstream and upstream of the genes. 
doi:l 0.1 371 /journa!.pone.0092502.t003 
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Figure 1. Distribution of minor allele frequencies (MAFs) of 
SNPs Identified from the T. rubirpes swimbladder. The X-axis 
represents the SNP minor allele frequency in percentage, while the Y- 
axis represents the number of SNPs with given minor allele frequency 
doi:1 0.1 371/journal.pone.0092502.g001 
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SNP distribution among genes and scaffolds 

SNPs distribution is important for consideration of coverage 
when using SNP makers. The distribution of SNPs in tlie genes 
was analyzed. Expressed short reads were mapped to a total of 
17,249 genes based on the fifth fugu T. rubripes genome assembly 
from Ensembl database. On average, 3.6 SNPs per gene were 
identified. A total of 1 1,306 expressed genes containing SNPs were 
identified in the swimbladder with the cutoff values of PRKM 
setting as 0.08. As shown in Figure 2, of these genes, 56.73% had 
fewer than 5 SNPs per gene. The number of genes with 26-30 
SNPs per gene is 40 and there are 30 genes harboring more than 
30 SNPs per gene. For instance, the dystonin (ENSTRUG 
00000015507) and annexin A5 (ENSTRUGOOOOOO 15464) have 
relatively large numbers of SNPs per gene, 73 and 63 SNPs, 
respectively. The fugu genome assembly (version 5.0) is composed 
of 7,1 19 scaffolds. The SNPs identified in the present study were 
found on the 2,246 scaffolds, about 27.7 SNPs per scaffold. As 
shown in Figure 3, a large number of scaffolds had fewer than 10 
SNPs per scaffold. The scaffold_l and scaffold_6 had the largest 
number of SNPs, 1,631 and 1,293 SNPs, respectively. 

Gene Ontology and KEGG pathway analysis 

Gene Ontology (GO) annotation was further performed for the 
annotated genes in terms of biological process, molecular function 
and cellular component. Distribution of the genes in different GO 
categories at level 2 is shown in Figure 4. In the swimbladder, 
8,922 expressed genes containing SNPs were assigned with one or 
more GO terms for biological process, molecular function and 
ceUular component. For biological process, genes involved in the 
metabolic process and cellular process were highly represented. 
For molecular function, binding was the most represented GO 
term, followed by catalytic activity. Regarding to the cellular 
component, the major categories were cell and cell part. The GO 
categories of expressed genes containing SNPs were in the same 
proportion to the GO categories of all the expressed genes 
(Figure 4). 

Besides GO analysis, KEGG pathway analysis was also carried 
out for the annotated genes, which is an alternative approach to 
categorize gene functions with the focus on biochemical pathways. 
A total of 3,808 expressed genes were assigned with one or more 
KEGG annotation and were mapped to KEGG pathways 
(Table 4). Of these annotated genes, 28.06% were classified into 
the Organismal Systems with the majority of which involved in 
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Figure 2. SNP distribution among genes. The X-axis represents 
gene size (number of SNPs per gene) 
doi:1 0.1 371/journal.pone.0092502.g002 
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Figure 3. SNP distribution among scaffolds. The X-axis represents 
scaffold size (number of SNPs per scaffold) 
doi:10.1371/journal.pone.0092502.g003 

immune system. Metabolism pathways including carbohydrate 
metabolism, amino acid metabolism and lipid metabolism 
represented 25.66%. Environmental information processing rep- 
resented 19.41%. The signal transduction was one of the weU- 
represented sub-pathways. In addition, 9.15% and 17.72% were 
classified into the Genetic information processing and Cellular 
Processes, respectively. 

Homologous genes containing SNPs between fugu 
swimbladder and human lung 

In this study, our KEGG pathway analysis indicated the tight 
junction existed, including 141 expressed genes containing SNPs. 
Tight junction is essential for epithelial morphology and 
function of swimbladder. Tight junctions serve to form seals 
among epithelial cells, creating a selectively permeable barrier to 
intercellular diffusion [27]. Claudins are transmembrane proteins 
which act in concert with other transmembrane and peripheral 
proteins to form the physical basis for tight junction [27,30]. In 
previous studies, claudin 4/5/6/7/9 genes were identified in the 
swimbladder of zebrafish [27] and 46 claudin genes in the fugu 
genome were identified and their phylogenetic relationships to 
those counterparts in mammals was determined [31]. In this study, 
16 members of claudin family were identified (Table 5). Three of 
the 16 claudin genes were highly expressed, including claudin 5a, 
5b and 7d. In the human airway, claudin 1, 3, 4, 5 and 7 are 
expressed in both bronchi and bronchioles. Claudin 5 is localized 
exclusively in the apical-most region of the tight junctions. Altered 
Claudin expression pattern can change the paracellular permeabil- 
ity characteristics of the epithelium. Claudin 5 overexpression 
increases the solute permeability [32,33]. Genome wide associa- 
tion studies showed the polymorphisms rs9290927, rs893051 and 
rsl 7501010 from clandin 1 were associated with nickel contact 
sensitization in individuals without ear piercings, contact sensiti- 
zation to fragrances, and with both organic compounds and nickel 
contact dermatitis in human, respectively [34]. The genetic 
variants in regulatory regions of clandin 1 can alter susceptibility 
to HCV infection [35]. 

In this study, 8 Wnt genes containing SNPs were identified and 
the expression levels of wnt 7b, wnt 5a and wnt 11 are higher 
(Table 6). Wnt signaling pathway has been reported to play 
important roles in mammalian lung development [36-38]. In 
previous studies, the down-regulation of Wnt signaling leading to 
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Figure 4. Gene Ontology of genes containing putative SNPs. 

doi:10.1371/journal.pone.0092502.g004 



defective swimbladder development in zebrafish was observed 
[39]. Wnt7h is expressed in the distal airway epithelium of lungs 
and plays critical roles in lung development such as distal epithelial 
cell fate decision, lung mesenchymal proliferation and smooth 
muscle differentiation [38,40-43]. It was found that wntSa is 
expressed in lung epithelium [38,44]. Wntll plays important roles 
in mouse lung development [38,45,46]. In chicken, 124 SNPs 
from 3 1 genes of Wnt signaling pathway were selected to genotype 
in 764 individuals resulted in 102 polymorphic SNPs [47]. In 
human, 14 SNPs from six Wnt pathway-related genes were 
genotyped in 210 individuals (145 men and 65 women), including 
Dickkopf 2 {DKK2) (rsl7037102, rs419558, and rs447372), DKK3 
(rs3206824, rsl 1022095, rsl472189, rs7396187, and rs2291599), 
DKK4 (rs2073664), sFRP4 (rsl802073 and rsl802074), SMAD7 
(rsl2953717), and DAAM2 (rs6937133 and rs2504106) [48]. Six 
common SNPs of WntlOb were identified in a sample of 1,029 
Korean female subjects, which were in almost complete linkage 
disequilibrium [49]. 

We observed the expression of two homologues of Ihh 
(ENSTRUGOOOOOO 12233 and ENSTRUGOOOOOO 13525) and 
Ptcl (ENSTRUGOOOOOO 145 14) containing SNPs from the swim- 
bladder transcriptome. The role of Hh (Hedgehog) signaling 
pathway in lung development is very crucial in human, mouse, 
chicken and Xenopus laevis [38,50-53]. Some development-related 
genes in lung had been identified in zebrafish, such as Sonic 
Hedgehog (^M)-related gene, Indian Hedgehog (/A/i)-related gene 
and their receptors, Patched l[Ptc T) and Ptc2 [54-61]. The human 
sonic hedgehog {SHH) gene is located in the 7q36 region, which is 
known to play an important role in embryo patterning, lung 
development and connection with sexual orientation. A SNP site 
(rs9333613) was found to be associated with male sexual 
orientation [62]. Ihh is a good candidate gene for association 
study of developmental disorders mainly affecting skeleton 
development. The previous study showed that the SNP sites of 



Ihh were associated with equine bone developmental disorders 
[63]. 

SNP validation 

As the SNPs reported in the present study were identified 
through bioinformatic analysis, the results were needed to evaluate 
for the validation rate. A total of 48 SNPs were randomly selected 
for validation by PGR amplification and Sanger sequencing [64] . 
Of the 48 SNPs, 26 SNPs (54%) were validated and 22 SNPs 
were not found by PGR amplification and direct sequencing 
(Table 7). 

Materials and Methods 

Ethics statement 

This study was approved by the Animal Care and Use 
committee of Key Laboratory of Mariculture & Stock Enhance- 
ment in North China's Sea at Dalian Ocean University. AU 
surgery was performed under sodium pentobarbital anesthesia, 
and all efforts were made to minimize suffering. 

Sample collection and RNA isolation 

A total of 45 Takifugu rubripes (length 20cm) were sampled from 
Dalian Tianzheng Industrial Co., Ltd (Dalian China). The 
swimbladders of these fish were collected and pooled. Tissues 
were placed into RNAlater (Ambion), stored at room temperature 
for 24 h, and then moved to -80°C for storage until RNA 
isolation. Total RNA was extracted from the pooled swim- 
bladder using the TRIzol R Reagent (Invitrogen, GA, USA) 
by following the manufacturer's protocol. The quantity and 
quality of total RNA was measured using an Agilent 2100 
Bioanalyzer. 
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Table 4. KEGG biochemical mappings for genes containing 
SNPs. 



Number of 

KEGG categories genes 
Metabolism 

Amino Acid Metabolism 344 
Biosynthesis of Polyl<etides and Nonribosomal Peptides 4 

Biosynthesis of Secondary Metabolites 69 

Carbohydrate Metabolism 422 

Energy Metabolism 162 

Glycan Biosynthesis and Metabolism 150 

Lipid Metabolism 331 

Metabolism of Cofactors and Vitamins 130 

Metabolism of Other Amino Acids 103 

Nucleotide Metabolism 239 

Xenobiotics Biodegradation and Metabolism 118 
Genetic information Processing 

Folding, Sorting and Degradation 266 

Replication and Repair 161 

Transcription 196 

Translation 116 
Environmental Information Processing 

Membrane Transport 39 

Signal Transduction 1142 

Signaling Molecules and Interaction 386 
Cellular Processes 

Behavior 28 

Cell Communication 522 

Cell Growth and Death 365 

Cell Motility 200 

Transport and Catabolism 316 
Organismal System 

Circulatory System 138 

Development 187 

Endocrine System 544 

Immune System 1116 

Nervous System 246 

Sensory System 35 

doi:1 0.1 371 /journal.pone.0092502.t004 



cDNA library construction and sequencing 

Total RNA was sent out for next generation sequencing 
provided by Beijing Institute of Genomics, Chinese Academy of 
Sciences. cDNA libraries were constructed from mRNA from 
swimbladder. cDNA libraries were prepared using the lUumina 
TruSeq RNA Sample Preparetion Kit (lUumina) according to the 
TruSeq protocol. After KAPA quantitation and dilution, the 
libraries were clustered 3 per lane and sequenced on an lUumina 
HiSeq 2000 instrument with 100 bp paired-end reads. 

Reads mapping 

The reads were mapping to the fugxi T. mbripes fifth genome 
assembly by BWA program. During the mapping phase, up to five 
mismatches were allowed. The expression levels (RPKM, Reads 
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Table 5. Identification of expressed Claudin genes containing 
SNPs. 



EnsembI Gene ID 


Gene name 


RPKM value 


Number of SNPs 


ENSTRUG00000018609 


Claudin 5b 


234.65 


5 


ENSTRUG00000016497 


Claudin 5a 


75.65 


5 


ENSTRUG00000007521 


Claudin 7a 


43.24 


6 


ENSTRUG00000010140 


Claudin 30c 


14.33 


3 


ENSTRUG00000004991 


Claudin 12 


8.98 


6 


ENSTRUG00000011829 


Claudin 1 la 


6.15 


4 


ENSTRUG00000015308 


Claudin 15a 


3.81 


1 


ENSTRUGOOOOaoa3031 


Claudin 25 


2.89 


2 


ENSTRUG00000010901 


Claudin 23 


1.83 


4 


ENSTRUG0000a001287 


Claudin 19 


1.63 


4 


ENSTRUG00000013204 


Claudin 18 


1.21 


2 


ENSTRUG00000007366 


Claudin 32a 


0.78 


3 


ENSTRUG00000011741 


Claudin 31 


0.76 


2 


ENSTRUG00000009832 


Claudin 28b 


0.31 


1 


ENSTRUG00000016459 


Claudin 15b 


0.21 


2 


ENSTRUG00000010378 


Claudin 23b 


0.14 


1 



doi:1 0.1 371 /journa!.pone.0092502.t005 



Per Kilobase of exon model per Million mapped reads) for each 
gene were calculated using uniquely mapped reads by in-house 
Perl script according to the equation: 

uniquci-eads X geneiengii, 

The cutoff value of gene expression was calculated for each 
sequencing sample, genes with RPKM greater than cutoff value 
were defined as expressed genes [65]. 

SNP identification 

BWA and SAMtools (Tools for alignments in the SAM format) 
software were used to align reads to the fugu genome assembly 
(version 5.0) for detecting SNPs [66,67]. Filtering thresholds were 



Table 6. Identification of expressed Wnt genes containing 
SNPs. 





EnsembI Gene ID 


Gene name 


RPKM value 


Number of SNPs 


ENSTRUG00000016453 


Wnt 7b 


48.59 


1 


ENSTRUG00000001530 


Wnt U 


29.03 


8 


ENSTRUG00000008614 


Wnt 5a 


28.52 


7 


ENSTRUG00000000172 


Wnt 4 


16.51 


4 


ENSTRUG00000003640 


Wnt 6 


11.5 


14 


ENSTRUG00000016522 


Wnt 5b 


1.71 


2 


ENSTRUG00000014284 


Wnt 9b 


1.38 


3 


ENSTRUG00000012568 


Wnt 2b 


0.33 


1 



doi:l 0.1 371 /journa!.pone.0092502.t006 
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Table 7. Summary of SNP validation. 





Enscmbl G6ne ID 


Gene Nsmfi 


Number of SNPs tested 


Number of SNPs vslidated 


ENSTRUGOOOOOOl 1 255 


Translocase of outer mitochondrial membrane 
20 homolog 


5 


2 


ENSTRUG00000008698 


RAB9A, member RAS oncogene family 


7 


4 


ENSTRUG00000014751 


Family with sequence similarity 46, member A 


5 


3 


ENSTRUG00000009192 


Coiled-coil domain containing 47 


5 


2 


ENSTRUG00000006299 


IViitochondrial ribosomal protein L21 


8 


5 


ENSTRUG00000006704 


Calpain small subunit 1 


5 


3 


ENSTRUG00000004026 


C-type lectin domain family 11, member A 


6 


3 


ENSTRUG00000014304 


Proliferating cell nuclear antigen 


7 


4 



doi:10.1371/journal.pone.0092502.t007 



set as: consensus quality is no less than 20 and coverage is no less 
than 10. 

Gene Ontology and KEGG pathway analysis 

Gene Ontology (GO) and KEGG pathway analyses were 
conducted to the genes containing SNPs. GO annotation analysis 
was performed using Blast2GO, an automated tool for the 
assignment of GO terms. The annotation result was categorized 
with respect to Biological Process, Molecular Function, and 
Cellular Component at level 2. In order to gain an overview of 
gene pathway networks, KEGG analysis was performed using the 
online KEGG Automatic Annotation Server (KAAS) (http:// 
www.genome.jp/kegg/kass/). The bi-directional best hit (BBH) 
method was used to obtain KEGG orthology assignments. 

SNP validation 

To evaluate the validation rate of the SNPs identified by 
bioinformatic analysis, we randomly selected 48 SNPs and 
validated by PGR amplification and direct sequencing. PGR 
primers were designed according to the assembled transcript 
sequences and were listed in the Table S2. Ten individuals were 
used for the SNP validation. 



Conclusions 

In this study, a large number of SNPs were identified by the 
transcriptome sequencing of the T. ruhirpes swimbladder using 
lUumina HiSeq2000 platform. A large proportion of randomly 
selected SNPs were verified using the Sanger sequencing, 
suggesting the high validation rate. The SNPs should provide 
valuable resources for genomic studies, evolution analysis, 
population genetic study, resource assessment, genetic hnkage 
analysis and genome-wide association studies. 

Supporting information 

Table SI The SNPs identified from the transcriptome of 
the swimbladder of Takifugu rubripes. 

(TXT) 

Table S2 Primers used for SNP validation in the study. 

(DOC) 
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